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ABSTRACT 


This  report  discusses  model  order  determination  algorithms  for  autoregressive 
filter  techniques  applied  to  tactical  radio  direction  finding.  In  particular  the  difficulties 
imposed  by  multipath  are  examined.  The  related  issue  of  the  determination  of  the  number 
of  signals  for  autoregressive  radio  direction  finding  techniques  using  the  eigen 
decomposition  method  is  also  discussed.  The  performance  of  both  model  order  and  signal 
number  determination  methods  are  also  evaluated  using  computer  simulations. 


RESUME 


Ce  rapport  contient  une  discussion  d’algorithmes  de  modeles  auto-regressifs  pour 
l’ordre  de  filtres  utilises  pour  la  determination  de  la  direction  d’arrivee  de  radios  tactiques. 
Les  difficultes  crees  par  les  trajectoires  multiples  recoivent  une  attention  speciale.  On 
discute  aussi  du  sujet  connexe  du  comptage  du  nombre  de  signaux  par  des  techniques  radio 
de  determination  de  direction  d’arrivee  basees  sur  des  vecteurs  propres.  Une  simulation  par 
ordinateur  est  utilisee  pour  evaluer  les  performances  des  methodes  de  determination  de 
l’ordre  et  du  nombre  de  signaux. 
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EXECUTIVE  SUMMARY 


This  report  reviews  methods  to  determine  the  optimum  model  order  of  adaptive 
filter  techniques  applied  to  the  problem  of  radio  direction  finding  (DF),  and  the  related 
issue  of  the  determination  of  the  number  of  signals  for  adaptive  filter  techniques  using 
eigen  decomposition  techniques.  This  work  was  carried  out  in  support  of  research  into 
advanced  radio  direction  finding  techniques  for  land  tactical  applications. 

For  radio  direction  finding  applications,  adjusting  the  adaptive  filter  model  order 
to  a  value  less  than  the  maximum  (i.e.  one  less  than  the  number  of  antennas  in  the  sensor 
array)  increases  the  ability  of  an  adaptive  filter  based  DF  algorithm  to  resolve  multipath 
signals.  It  is  from  this  point  of  view  that  model  order  determination  methods  are  discussed 
in  this  report. 

In  general,  the  optimum  model  order  is  not  known  beforehand,  so  a  number  of 
model  order  determination  methods  have  been  proposed  to  tackle  this  problem.  In  this 
report,  these  are  classified  according  to  the  autocorrelation  matrix  estimator  used.  Three 
classes  of  estimators  are  considered,  which  are:  (a)  the  autocorrelation  method,  (b)  the 
covariance  and  modified  covariance  method,  and  (c)  the  eigen  decomposition  method. 

The  autocorrelation  method  is  not  appropriate  for  accurate  direction  finding 
purposes  since  adaptive  filter  DF  techniques  based  on  this  method  have  very  poor 
resolution.  However,  it  is  discussed  in  this  report  since  the  first  model  order  determination 
algorithms  were  proposed  for  this  method,  and  many  of  the  problems  encountered  are 
common  to  all  classes.  The  performance  of  these  algorithms  is  quite  poor,  typically  picking 
model  orders  that  are  too  low. 

The  covariance  and  modified  covariance  methods  are  much  superior  for  direction 
finding  purposes  than  the  autocorrelation  method.  It  has  been  determined  by  researchers 
[10]  both  analytically  and  experimentally  that  a  model  order  of  between  Nj 3  and  N/2  is 
suitable,  where  N  is  the  number  of  antenna  sensors.  Model  orders  higher  than  this  lead  to 
spurious  estimates  which  severely  degrade  bearing  accuracy. 

The  eigen  decomposition  method  enhances  the  signal  to  noise  ratio  of  the 
autocorrelation  matrix  estimated  using  either  the  covariance  or  modified  covariance 
method.  This  enhancement  suppresses  spurious  bearing  estimates  so  that  higher  model 
orders  can  be  used.  Values  of  between  27V/3  and  3A/4  have  been  proposed. 

An  additional  parameter  required  to  optimally  enhance  the  signal  to  noise  ratio 
using  the  eigen  decomposition  method  is  the  number  of  radio  signals  present.  Since  this 
problem  is  closely  related  to  the  model  order  determination  problem,  model  determination 
methods  have  been  adapted  for  this  purpose.  For  widely  spaced  (in  bearing)  signals,  and 
signal  to  noise  ratios  of  20  dB,  the  adapted  algorithms  work  extremely  well.  Under  more 
realistic  conditons  (i.e.  closely  spaced  signals,  lower  signal  to  noise  ratios,  etc.)  it  is 
expected  that  the  ability  of  these  signal  number  determination  algorithms  to  accurately 
determine  the  number  of  signals  will  degrade,  but  not  before  the  bearing  estimation 
accuracy  of  DF  algorithms  using  the  eigen  decomposition  method  also  degrade. 

In  general,  methods  to  determine  the  number  of  signals  used  in  conjunction  with 
adaptive  filter  DF  methods,  where  both  are  based  on  using  the  eigen  decomposition 
method,  appear  to  offer  the  best  performance.  The  worst  performance  is  provided  by  model 
order  determination  methods  used  in  conjunction  with  DF  methods  where  both  are  based 
on  the  autocorrelation  method.  In  this  case  performance  is  very  poor,  and  inappropriate  for 
direction  finding  purposes  when  accurate  bearing  estimation  is  necessary. 
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1.0  INTRODUCTION 


Current  conventional  land  tactical  radio  direction  finding  systems  suffer  from  two 
major  problems.  One  is  the  inability  to  accurately  determine  the  bearing  of  several  radio 
transmitters  simultaneously,  and  the  second  is  a  degradation  of  performance  in  the 
presence  of  multipath.  A  considerable  amount  of  research  and  development  has  been 
carried  out  over  the  last  twenty  years  in  advanced  spectral  estimation  algorithms  which 
has  led  directly  to  the  development  of  superresolution  radio  direction  finding  (DF) 
algorithms  which  may  be  capable  of  dealing  with  these  problems.  In  particular,  adaptive 
filter  methods  (including  eigen  based  methods)  have  shown  promise  in  this  area  with 
performances  approaching  maximum  likelihood  methods  [1]. 

The  superresolution  abilities  of  adaptive  filter  methods  is  based  on  matching  the 
sensor  data  with  an  appropriate  signal  model.  At  the  heart  of  these  methods  lies  a  system 
of  equations,  formed  from  the  antenna  sensor  data,  whose  solution  provides  the  signal 
bearings.  If  M  signals  impinge  upon  the  sensor  array  (where  Mis  less  than  the  number  of 
sensors),  then  M  linearly  independent  equations  are  required  to  solve  for  the  M  signal 
bearings. 


In  the  case  where  the  signals  are  all  noncoherent  (transmitted  by  independent 
transmitters)  M  linearly  independent  equations  may  be  formed  from  M  or  more  sensor  data 
vectors.  In  this  report,  a  data  vector  is  defined  as  a  single  set  of  measurements  of  all  the 
sensor  outputs  taken  at  time  instant  tm.  Data  vectors  are  assumed  to  be  u.  formly  sampled 
in  time  where  the  sampling  rate  is  low  enough  so  that  consecutive  data  vectors  are  not  fully 
correlated  (i.e.  not  scalar  versions  of  each  other). 

In  the  case  of  multipath  signals,  consecutive  data  vectors  remain  correlated  in 
time  no  matter  what  the  sampling  interval  (except  for  the  temporal  noise  component)  so 
that  effectively  only  one  data  vector  is  available.  To  generate  the  required  number  of 
linearly  independent  equations,  the  data  is  divided  into  subsets  or  subarrays  from  which 
the  required  number  of  linearly  independent  equations  may  be  generated.  This  is  discussed 
in  more  detail  in  section  1.2.  In  adaptive  filtering  terminology,  the  number  of  antenna 
sensors  in  the  subarray,  minus  one,  is  called  the  model  order.  In  direction  finding 
terminology,  the  technique  of  using  lower  model  orders  to  generate  more  linearly 
independent  equations  is  called  spatial  averaging  [2]. 

From  the  preceding  discussion,  it  is  apparent  that  spatial  averaging  is  a  necessity 
when  direction  finding  on  multipath  signals,  and  completely  unnecessary  if  the  signals  are 
noncoherent  and  a  sufficient  number  of  data  vectors  (>M)  have  been  measured. 
Consequently,  model  order  selection  is  discussed  in  the  context  of  determining  the  bearing 
of  the  direct  and  indirect  paths  of  a  signal  in  a  multipath  environment. 

A  number  of  methods  to  determine  model  order  and  the  number  of  signals  present 
are  evaluated  in  this  report.  In  particular,  they  are  evaluated  with  respect  to  the 
limitations  imposed  by  multipath  and  tactical  constraints  which  restrict  the  number  of 
sensors.  The  effect  of  multipath  is  equivalent  to  DF  processing  using  a  single  data  vector 
since  additional  vectors  are  fully  correlated  (ignoring  the  noise).  Consequently  only  single 
vector  processing  is  considered  in  this  report.  To  satisfy  the  tactical  constraint  the  sensor 
arrays  are  limited  to  16  antenna  sensors.  Additionally,  since  spatial  averaging  can  only  be 
used  for  linear  sensor  arrays  with  uniformly  spaced  sensors,  only  these  type  of  sensor  arrays 
will  be  considered. 
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Noise  discussed  in  this  report  is  assumed  to  be  additive  white  Gaussian  in  nature 
and  uncorrelated  between  sensors.  For  a  single  data  vector  the  effects  of  temporal  and 
spatial  noise  are  indistinguishable  so  no  distinction  is  made  between  these  two  types  of 
noise. 

1.1  Spectral  Estimation  and  Radio  Direction  Finding 

The  application  of  spectral  estimation  concepts  to  radio  direction  finding  is 
illustrated  using  Figure  1.  A  single  radio  signal  in  a  noiseless  environment  impinges  on  an 
antenna  baseline  with  some  arbitrary  direction  of  arrival.  The  amplitude  and  phase 
measured  at  sensor  n  is  represented  by  the  complex  value, 

-;(nus  d>8) 

2n  =  y[t)e  (1.1) 

where  y(i)  represents  the  complex  baseband  modulation  envelope  (for  all  types  of 
modulation  and  is  assumed  to  remain  constant  for  the  measurement  of  a  single  vector),  d  is 
the  spacing  between  consecutive  sensors,  0\s  the  signal  phase  which  is  a  function  of  the 
distance  from  the  transmitter  to  sensor  0,  and 

u)s  =  -j-  cos(0)  (1-2) 


Here  A  is  the  wavelength,  and  (f>  is  the  signal  bearing  in  azimuth.  (Signals  with  elevation 
angles  different  than  0  degrees  are  not  considered  here).  Inspection  ol  equation  (1.1)  reveals 
that  in  represents  vectors  of  a  complex  sinusoid  measured  spatially  (i.e.  with  respect  to  the 
sensor  baseline)  and  that  u,  is  the  spatial  frequency  of  the  complex  sinusoid. 


FIGURE  1:  Radio  direction  finding  sensor  array 


If  a  number  of  signals  with  different  bearings  impinge  upon  the  sensor  system  then 
the  spatial  signal  xn  is  the  sum  of  a  number  of  complex  sinusoids  of  different  amplitudes 
and  frequencies  of  the  form  of  equation  (1.1).  Spectral  estimation  techniques  can  be  used  to 
resolve  the  various  components  by  first  generating  the  spatial  power  spectral  density 
function.  By  taking  advantage  of  equation  (1.2),  this  spectrum  is  converted  to  a 
DF  spectrum  which  gives  the  power  density  versus  bearing.  The  location  of  the  peaks  (local 
maximums)  in  the  DF  spectrum  are  then  used  as  the  estimates  for  the  signal  bearings. 
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1.2  Autoregressive  Model 

The  discussion  that  follows  focuses  on  adaptive  filtering  spectral  estimation 
techniques  which  have  been  adapted  for  the  purposes  of  radio  direction  finding.  The  types 
of  filters  used  in  these  applications  can  be  classified  as  autoregressive  (ART  moving 
average  (MA),  or  autoregressive-moving  average  (ARMA).  In  general,  techniques  involving 
either  MA  or  ARMA  filters  are  computationally  more  difficult  since  they  involve  a  set  of 
nonlinear  equations,  and  for  the  purposes  of  this  report  will  not  be  considered. 

Autoregressive  methods  are  based  on  the  idea  that  a  uniformly  sampled  sequence 
of  complex  antenna  sensor  data,  {xo,  ii,  X2,  ...,  x^-i},  (where  iV  represents  the  number  of 
sensors)  can  be  modelled  as  a  linear  rational  system  driven  by  white  noise. 

Mathematically, 


%m-n  “I"  % 


mt 


(1.3) 


where  v,  is  a  white  noise  source  with  a  variance  of  a2,  and  p  is  the  model  order  which  is 
chosen  to  be  less  than  N.  For  a  particular  model  order  p,  the  values  of  the  tap  weight 
coefficients  a„  are  usually  chosen  so  that  the  variance  of  the  white  noise  vt  is  minimized. 
Equation  (1.3)  also  represents  a  linear  predictor  where  p  past  data  vectors  are  used  to 
predict  the  current  value  xm.  In  this  context,  vm  represents  the  prediction  error. 


Based  on  minimizing  the  noise  variance,  or  prediction  error,  the  optimum 
weighting  coefficients  can  be  found  by  solving  the  augmented  normal  matrix  equations 
given  by 


Ra  =  a2u, 

where  the  (p+l)*(p+l)  autocorrelation  matrix  Ris  defined  by, 

R  = 

♦ 

'  E{x o  xo}, 
E{x i  xo}, 
E{x2  x0}, 

*  * 

E{x o  xi},  E{x o  12},  .. 

E{xi  xi},  E{x  1  x2} ,  .. 

E\X  2X1},  E{x  2X2},.. 

* 

.,  E{x 0  x^}  ' 
.,  £{x,  x^} 

.,  E{x 2  ip} 

•  * 

.  E{xpx0}, 

•  *  •  * 
E{xpxi}  E{xp  x2}  .. 

•  * 

. ,  E{  xp  xp}  _ 

the  augmented  tap  weight  vector  a  is  defined  by, 


(1.4) 


(1.5) 


1 

Ol 


(1.6) 


3 


and  the  vector  u  is  defined  by, 


‘  1 ' 
0 


0 


(1.7) 


Using  autoregressive  filter  techniques,  an  autoregressive  DF  estimator  can  be 
defined  by  [3] 

S(<f>)  =  7^2.  (1-8) 

|  aHe  | 

where  e  is  the  p+1  element  steering  vector  defined  by, 


e  = 


1 

*io), 2d 

e  a 


(1.9) 


e*jo>spd 


and  us  was  defined  previously  by  equation  (1.2). 
1.3  The  Estimated  Autocorrelation  Matrix 


To  this  point,  the  calculation  of  the  tap  weight  filter  coefficients,  and  the 
formation  of  the  spatial  spectral  density  function  (equation  (1.8))  using  these  coefficients 
has  been  discussed.  The  assumption  was  made  that  the  optimum  model  order  p  was 
known.  In  the  case  where  the  true  autocorrelation  matrix  is  known,  the  higher  the  order  of 
p,  the  better  the  resolution  of  the  AR  spectral  estimator.  Consequently,  when  the  true 
autocorrelation  matrix  is  known,  the  best  value  of  p  is  equal  to  N- 1. 

In  most  practical  estimation  problems  the  true  autocorrelation  matrix  is  unknown 
and  must  be  estimated  from  the  data.  The  estimated  autocorrelation  matrix  can  be 
defined  by  the  equation 

R  =  £x,x»  +  ^yj7 )  (1.10) 

* 

where  x*  represents  a  p+1  length  column  vector  containing  the  sensor  data  in  the  forwards 
direction.  That  is, 


and  y  represents  a  p+1  length  column  vector  containing  the  conjugated  sensor  data  in  the 
reverse  direction.  That  is, 


* 


y  i  = 


xi*p 
xi*p- 1 
xi\v 2 


L  xi  J 


(1.12) 


The  range  of  the  indices  i  and  j  have  purposely  been  left  unspecified  since  different 
methods  use  different  ranges. 


2.0  MODEL  ORDER  DETERMINATION  FOR  THE 
AUTOCORRELATION  METHOD 

The  autocorrelation  method  of  estimating  the  autocorrelation  matrix  can  be 
defined  from  equation  (1.10)  as 


-  „ 

R  =  ]>\xt  (2.1) 

»*-  p 

Since  this  expression  involves  unknown  data  (z.i,  x.2,  etc.,  and  xN,  Xn*  1,  etc.),  the  unknown 
data  is  assumed  to  be  zero.  In  other  words,  the  data  is  pre-  and  postwindowed.  Windowing 
results  in  a  loss  of  resolution  which  becomes  progressively  worse  as  the  number  of  sensors, 

N,  decreases.  As  a  result,  the  autocorrelation  method  is  generally  not  favoured  for  spectral 
estimation  problems  involving  short  data  records  (such  as  direction  finding),  however,  it  is 
mentioned  here  since  some  of  the  earliest  model  order  determination  algorithms  were 
evaluated  using  this  approach. 

Figure  2  is  an  example  of  the  effect  that  model  order  selection  has  on  the  AR 
spectral  estimate  where  the  autocorrelation  matrix  described  by  equation  (1.10)  was  used. 
In  this  example  a  radio  direction  finding  system  using  a  linear  16  element  sensor  array  with 

O. 5A  spacing  was  simulated.  The  actual  signals  had  unit  power  and  bearings  of  40  and  100 
degrees  respectively.  The  noise  variance  was  0.01. 

Inspection  of  Figure  2  reveals  that  when  the  model  order  is  less  than  the  number 
of  signals  (in  this  case  p<2),  the  resolution  of  the  spectral  estimator  is  seriously  degraded. 
In  this  case  there  are  fewer  filter  coefficients  than  there  are  signals. 

When  the  model  order  becomes  too  large  (in  this  case  p>4)  the  spectral  peaks 
corresponding  to  the  signal  locations  become  corrupted  and  spurious  peaks  begin  to  appear. 
In  this  case  there  are  more  filter  coefficients  than  there  are  signals.  This  leads  to  better 
resolution,  however,  the  extra  filter  coefficients,  which  model  the  noise,  occasionally  give 
rise  to  spurious  estimates  as  seen  in  this  example. 

Another  problem  is  that  for  successively  higher  model  orders  the  elements  of  the 
autocorrelation  matrix  are  formed  from  fewer  and  fewer  data  values.  This  leads  to  greater 
variance  in  the  autocorrelation  matrix  estimate  and  in  turn  leads  to  greater  variance  in  the 
spectral  estimate. 
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(i)  Modal  Order  9 
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FIGURE  2:  Effect  of  the  choice  of  AR  model  order  on  the  estimated  DF  spectrum. 

Actual  signal  bearings  are  shown  by  the  dashed  lines.  Signal  powers 
(|  amplitude  I2)  were  1.0  and  the  noise  variance  was  0.01  (i.e.  signal  to  noise 
ratio  of  20  dB). 
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To  quantify  the  effect  of  model  order  selection  on  DF  accuracy,  100  DF  trials  were 
simulated  for  a  two  signal  case  (M=  2)  and  a  five  signal  case  (M=  5)  using  the  same  DF 
sensor  system  described  for  Figure  2  (i.e.  16  elements).  The  results  are  shown  in  Figure  3. 
For  each  trial  the  signal  powers  were  set  to  1,  the  signal  phases  were  randomly  chosen 
between  0  and  2n,  and  the  signal  bearings  were  set  to  40  and  100  degrees  for  the  two  signal 
case,  and  40,  70,  100,  140,  and  160  for  the  five  signal  case.  The  noise  variance  was  set  at 
0.01  (SNR  =  20  dB).  Bearing  errors  were  calculated  by  locating  the  M  largest  peaks  in  the 
DF  spectrum  and  subtracting  these  bearings  from  the  true  bearings.  In  cases  where  there 
were  less  than  M  peaks  in  the  DF  spectrum,  such  as  when  the  model  order  p  <  M,  the  error 
was  calculated  as  the  difference  between  the  corresponding  estimated  bearing  and  the 
nearest  true  bearing.  The  fact  that  some  bearings  were  completely  missed  this  way  was  not 
taken  into  account  in  the  results.  The  accuracy  measurements  shown  in  Figure  3  were 
computed  as  the  standard  deviation  of  the  bearing  error  for  a  particular  model  order. 


Bearing  Accuracy 
(degrees) 


(b) 


Model  Order  (p) 


FIGURE  3:  Bearing  estimation  accuracy  versus  model  order  for  the  case  of 

equi-amplitude  signals  (SNR  =  20  dB)  with  bearings  of:  (a)  40  and  100 
degrees,  and  (b)  40,  70,  100,  140,  and  160  degrees. 


From  the  results  shown  in  Figure  3,  it  is  clear  that  the  optimum  model  order  is 
equal  to  or  slightly  greater  than  the  number  of  signals  present.  The  large  error  that  occurs 
for  model  order  selections  that  are  too  low  or  too  high  supports  the  observations  made  in 
Figure  2.  It  should  be  pointed  out  that  Figure  3  is  highly  contrived  (compared  to  what 
might  be  expected  in  a  true  multipath  environment)  since  the  parameters  of  the  true 
signals  used  in  Figure  3  were  chosen  to  be  within  the  resolution  capabilities  of  the  AR 
spectral  estimator,  i.e.  bearings  were  widely  spaced,  and  signal  powers  were  reasonably 
strong  compared  to  the  noise  power.  However,  in  cases  where  the  AR  estimator  is  unable  to 
resolve  two  or  more  signals,  the  accuracy  of  the  estimate  degrades  rapidly  to  the  point  that 
the  estimates  for  any  model  order  are  not  much  more  accurate  than  random  guesses.  In 
general,  Figures  2  and  3  indicate  the  need  for  some  method  to  choose  the  correct  model 
order. 


Algorithms  designed  to  determine  the  best  model  order  take  advantage  of  two 
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properties  of  AR  filters.  The  first  property  is  that  as  the  model  order  is  increased  from  1  to 
N- 1,  the  rate  of  decrease  of  the  prediction  error  variance  a 2  (see  equation  (1.4)),  is  greater 
when  the  model  order  is  less  than  the  number  of  signals  than  compared  to  when  the  model 
order  is  greater  than  the  number  of  signals.  The  second  property  is  that  as  the  model  order 
increases,  the  variance  in  the  estimates  of  the  AR  tap  weight  coefficients  increases. 

30 
25 
20 
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Variance  in  dB  10 
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FIGURE  4:  Prediction  error  variance  versus  model  order  for  the  example  shown  in 
Figure  2. 


Figure  4  demonstrates  the  first  property  where  the  "kink"  in  the  prediction  error 
variance  curve  corresponds  to  the  number  of  signals,  which  in  this  case  was  2.  Although 
one  might  be  tempted  to  base  model  order  determination  on  locating  this  kink,  the  problem 
becomes  far  more  difficult  when  closely  spaced  signals  in  bearing,  lower  signal  to  noise 
ratios  and  signals  of  different  powers  are  considered,  as  illustrated  in  Figure  5.  However,  by 
taking  advantage  of  the  second  property,  the  optimum  model  order  can  be  estimated  by 
comparing  the  decreasing  prediction  error  (as  shown  in  Figure  4)  versus  increasing  bearing 
estimation  variance  as  the  model  order  is  increased.  Most  of  the  model  order  determination 
algorithms  discussed  in  this  report  are  based  on  this  idea. 
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FIGURE  5:  Prediction  error  variance  versus  model  order  for  signal  bearings  of  40  and  70 
degrees,  and  powers  of  1.0  and  0.1  units. 
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One  of  the  earliest  approaches,  called  the  Final  Prediction  Error  algorithm  (FPE) 
[4]  was  proposed  by  Akaike  in  1969.  It  is  given  by, 


FPE(p)  =  a) 


(2.2) 


where  is  the  variance  of  the  prediction  error  for  the  AR  model  order  p.  It  is  assumed 
that  the  sample  mean  has  been  subtracted  from  the  data  (this  assumption  also  applies  to 
the  rest  of  the  algorithms  discussed  in  section  2.0).  The  minimum  of  the  above  function 
gives  the  estimate  of  the  optimum  model  order. 

In  1974  Akaike  proposed  another  model  order  selection  algorithm  called  the 
Akaike  Information  Criterion  (AIC)  [4]  based  on  a  maximum  likelihood  approach. 
Assuming  the  process  has  Gaussian  statistics,  the  algorithm  is  given  by 


AIC(p)  =  Nl*(o§  +  2p  (2.3) 

The  model  order  is  again  selected  as  the  value  which  minimizes  this  function. 

Statistically,  it  was  determined  that  the  AIC  method  was  inconsistent  since  the 
probability  of  error  in  choosing  the  correct  model  order  does  not  approach  zero  as  the 
number  of  data  values  A  approaches  infinity.  In  response  to  this,  Rissanen  developed  a 
variant  algorithm  in  1983  which  is  statistically  consistent  for  large  N,  called  the  Minimum 
Description  Length  (MDL)  [5].  Again  assuming  Gaussian  statistics,  the  model  order  is 
determined  by  minimizing  the  following  function, 


MDL(p)  =  JVln(ap)  +  2pln(A)  (2.4) 

At  the  same  time  that  Akaike  introduced  the  AIC  algorithm,  Parzen  introduced 
an  algorithm  he  called  the  Criterion  Autoregressive  Transfer  function  (CAT)  [6].  The  best 
estimate  of  the  optimum  model  order  is  determined  by  minimizing  the  function 


CAT(p)  -  (-JJ  ^  J';1  )  +  Vpl 
;  =  i 


(2.5) 


where 


ui  =  [N/{N-  j)\  <7;  (2.6) 

Figure  6  provides  two  examples  of  the  performance  of  these  algorithms  for  a  single 
sensor  data  vector  (again  using  a  16  element  array  with  the  same  signal  and  noise 
parameters  used  in  Figure  3).  In  these  two  examples,  all  the  algorithms  were  successful  in 
choosing  the  optimum  model  order  (i.e.  equal  to,  or  slightly  greater  than  the  number  of 
signals). 
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FIGURE  6:  Model  order  determination  algorithms  for  a  (a)  two  signal  case,  and  (b)  a 
five  signal  case 


In  general,  the  derivation  of  the  algorithms  described  in  this  section  depend  on 
approximations  which  are  true  for  large  samples  (i.e.  N  large)  but  do  not  necessarily  hold 
in  the  case  of  short  data  records.  Researchers  have  found,  that  for  the  case  of  short  data 
records  involving  sinusoids  in  noise,  none  of  the  algorithms  work  well  [7]  and  in  general 
tend  to  pick  a  model  order  that  is  too  low.  This  conclusion  is  also  supported  by  the  results 
shown  Figure  7  which  gives  a  histogram  of  the  predicted  model  orders  for  each  algorithm 
for  100  trials  of  a  5  signal  environment  (the  same  sensor  data  that  was  used  in  Figure 
3(b)). 


10 


Model  Order 


FIGURE  7:  Performance  of  the  model  order  determination  algorithms  using  the  sensor 
data  from  Figure  3(b),  i.e.  100  trials  of  a  16  element  system  in  a  five  signal 
environment. 
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3.0  MODEL  ORDER  DETERMINATION  FOR  THE  COVARIANCE 
AND  MODIFIED  COVARIANCE  METHOD 

The  covariance  method  of  estimating  the  autocorrelation  matrix  uses  only  known 
data,  thus  avoiding  the  use  of  windowing  which  degrades  the  resolution  of  the  spectral 
estimator.  Based  on  equation  (1.10)  the  covariance  method  can  be  defined  as, 

.  H 

R=  ^xpc,.  (3.1) 

i  =o 

An  important  feature  of  the  autocorrelation  matrix  estimated  in  this  manner  is  that  the 
rank  of  the  matrix  is  the  lesser  of  N-p  (the  number  of  linear  independent  data  vectors 
used),  or,  p+ 1,  (the  length  of  the  data  vector). 

The  modified  covariance  method  takes  advantage  of  the  relationship  between  the 
forward  data  vector  x<  and  the  conjugated  reverse  data  vector  y,  given  by, 

£{*,  4  =  £{y,  y!}.  (3.2) 

This  is  used  to  generate  an  improved  autocorrelation  matrix  estimate  which  can  be  defined 
as 


R=  +  Ay,yi  (3-3) 

«  *0  t  =0 

The  advantage  of  this  method,  compared  to  the  covariance  method,  is  that  each 
autocorrelation  matrix  element  is  produced  from  twice  as  many  data  values  which  results 
in  less  variance  in  the  estimates.  Additionally  since  twice  as  many  linearly  independent 
data  vectors  are  used,  the  rank  of  the  matrix  is  the  lesser  of  2 (N-p)  or  p+1.  Based  on  rank 
considerations,  the  maximum  number  of  signal  bearings  that  may  be  solved  using  the 
modified  covariance  method  is  2N/3  as  compared  to  N/2  for  the  covariance  method.  Given 
these  advantages,  and  since  very  little  extra  processing  is  involved,  the  modified  covariance 
method  is  the  more  preferred  method  for  direction  finding  type  problems.  Most  of  the 
following  discussion  centers  on  the  modified  covariance  method.  However,  the  remarks 
made  concerning  this  method  may  readily  be  extended  to  the  covariance  method  by  taking 
the  differences  discussed  above  into  account. 

Figure  8  illustrates  the  effect  of  changing  the  model  order  on  the  spectral 
estimates  when  the  modified  covariance  method  is  used  (the  sensor  data  in  this  example 
was  the  same  used  in  Figure  2).  As  in  the  autocorrelation  method,  the  best  results  occur  for 
model  orders  equal  to  or  slightly  greater  than  the  number  of  signals  present. 
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FIGURE  8:  Effect  of  the  choice  of  AR  model  order  on  the  estimated  DF  spectrum. 

Actual  signal  bearings  are  shown  by  the  dashed  lines.  Signal  powers  were  1.0 
and  the  noise  variance  was  0.01  (i.e.  signal  to  noise  ratio  of  20  dB). 
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Figure  9  shows  the  effect  of  changing  model  order  on  the  accuracy  of  the  bearing 
estimates  using  the  modified  covariance  method.  The  sensor  data,  and  method  of 
computing  the  sensor  data  was  described  previously  in  section  2.0  for  Figure  3.  An 
examination  of  Figure  9(b)  shows  that  this  plot  can  be  broken  down  into  5  separate 
regions.  In  the  first  region,  p<5,  the  model  order  is  less  than  the  number  of  signals  present 
so  that  it  is  not  possible  to  determine  all  the  signal  bearings.  As  in  the  autocorrelation 
method,  this  results  in  poor  estimates  of  the  spatial  frequency  spectrum. 


Model  Order  (p) 


FIGURE  9:  Bearing  estimation  accuracy  versus  model  order  for  the  case  of 

equi-amplitude  signals  (SNR  =  20  dB)  with  bearings  of:  (a)  40  and  100 
degrees,  and  (b)  40,  70,  100,  140,  and  160  degrees. 


In  the  second  region,  5<p<7,  the  model  order  is  sufficiently  high  that  it  is  now 
possible  to  determine  all  the  signal"  bearings,  but  not  high  enough  that  spurious  estimates 
occur.  In  comparing  the  estimation  accuracy  in  this  region  to  the  rest  of  the  plot,  it  is  clear 
that  this  region  represents  the  optimum  order  for  the  AR  spectral  estimator. 

In  the  third  region,  8<p<  10,  spurious  estimates  occur  which  degrade  the  accuracy 
of  the  bearing  estimates  (i.e.  if  a  spurious  peak  in  the  DF  spectrum  is  larger  than  any  of 
the  signal  peaks,  it  is  then  mistakenly  assumed  to  be  a  signal  peak).  The  error  introduced 
by  the  spurious  estimates  increases  as  the  number  of  extra  model  coefficients  ( p-M) 
increases.  Additionally  as  the  order  increases,  the  variance  in  the  estimates  increases  as 
well  due  to  the  decreasing  number  of  data  values  used  to  form  each  element  of  the 
autocorrelation  matrix. 

In  the  fourth  region,  ll<p<13,  although  the  model  order  is  increasing,  the  rank  of 
the  autocorrelation  matrix  decreases  according  to  the  function  2 (N-p).  As  a  result,  the  filter 
tap  weight  coefficients  can  no  longer  be  uniquely  determined.  An  appropriate  choice  for  the 
tap  weight  coefficient  vector  (defined  by  equation  (1.6))  is  the  vector  which  has  the 
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minimum  norm  (where  the  vector  norm  is  given  by  aHa)  of  all  possible  solutions.  This 
vector  is  unique  and  has  the  desirable  property  that  spurious  estimates  are  suppressed  [8]. 
Using  the  tap  weight  coefficients  generated  in  this  manner,  the  reduced  rank  effectively 
reduces  the  number  of  signals  (including  spurious  signals)  that  may  be  uniquely 
determined.  In  consequence,  as  the  order  increases,  the  rank  decreases,  and  the  number  of 
spurious  signals  and  resultant  bearing  error  variance  decreases. 

A  special  case  occurs  when  the  rank  of  the  matrix  given  by  2 (N-p),  equals  the 
number  of  signals  A/(such  as  for  p=15  in  Figure  9(a)).  This  is  called  tne  Kumaresan-Prony 
case  [9]  and  is  of  special  interest  since  only  the  actual  signal  bearings  can  be  uniquely 
determined  i.e.  spurious  signals  do  not  occur.  Choosing  the  model  order  corresponding  to 
this  special  case  is  also  equivalent  to  using  noise  suppression  techniques  discussed  in  section 
4.0  without  the  penalty  of  additional  processing.  The  difficulty  is  that  knowledge  of  the 
exact  number  of  signals  is  required,  an  issue  which  is  also  discussed  in  section  4.0. 

In  the  fifth  region,  p>13,  the  rank  of  the  autocorrelation  matrix  decreases  below 
the  minimum  required  value  (in  this  case  5)  and  so  suffers  the  same  problems  as  those  of 
the  first  region,  plus  in  addition,  the  effect  of  increasing  variance  in  the  autocorrelation 
matrix  estimate. 

To  determine  the  correct  model  order,  the  algorithms  discussed  in  section  2.0 
could  be  used.  However,  as  discussed  previously,  these  algorithms  have  not  been  found  to 
work  well  for  short  data  segments.  Lang  and  McClellan  [10]  showed  analytically  that  a 
model  order  selection  of  about  N/ 3  provides  good  results  for  both  the  covariance  and 
modified  covariance  methods.  This  seems  reasonable  in  view  of  the  results  shown  in 
Figure  9,  assuming  the  number  of  signals  M<N/3,  otherwise  errors  will  occur  since  there 
would  be  an  insufficient  number  of  tap  weight  coefficients  to  estimate  all  the  signal 
bearings. 


4.0  MODEL  ORDER  AND  SIGNAL  NUMBER  DETERMINATION  USING 
THE  EIGEN  DECOMPOSITION  METHOD 

Improved  AR  spectral  estimates  can  be  made  by  taking  advantage  of  eigen 
decomposition  techniques  to  generate  a  better  estimate  of  the  noiseless  autocorrelation 
matrix,  and  improve  the  spectral  estimates.  A  summary  of  this  method  is  given  in  the 
following  paragraphs. 

For  the  case  of  signals  in  white  noise,  the  autocorrelation  matrix  can  be 
represented  as 


M-  1 

R=^A*vm  +  ^2I,  (4.1) 

*  =  o 

where  A*  represents  the  nonzero  eigenvalues  of  the  noiseless  autocorrelation  matrix  and  are 
ordered  from  largest  (Ao)  to  smallest,  v*  represents  the  corresponding  eigenvectors  and  rj 2 
represents  the  noise  variance  of  the  sensor  noise.  The  first  term  represents  the  signal 
correlation  matrix  where  it  is  assumed  that  M<p  and  the  second  term  represents  the  noise 
correlation  matrix.  Ideally,  knowing  the  signal  correlation  matrix  would  allow  the 
calculation  of  the  exact  signal  bearings. 
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FIGURE  10:  Effect  of  the  choice  of  AR  model  order  on  the  estimated  DF  spectrum. 

Actual  signal  bearings  are  shown  by  the  dashed  lines.  Signal  powers  were  1.0 
and  the  noise  variance  was  0.01  (i.e.  signal  to  noise  ratio  of  20  dB). 
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The  actual  eigen  decomposition  of  the  autocorrelation  matrix  is  given  by, 

w- 1  N-l 

R  =  ^(A*+?72)v*v!  +  ^Vv*v2  (4.2) 

*= 0  k-U 

where  A*+f?2  represents  the  M largest  eigenvalues  of  R  and  r?2  represents  the  N-M+l 
smallest  eigenvalues  of  R.  The  first  summation  term  in  the  above  expression  is  called  the 
signal  subspace  since  it  is  spanned  by  the  signal  vectors  (plus  noise)  while  the  second 
summation  term  is  called  the  noise  subspace  since  it  is  spanned  only  by  noise  vectors. 

The  above  results  are  only  approximately  true  in  the  case  of  the  estimated 
autocorrelation  matrix,  and  only  apply  if  either  the  covariance  or  modified  covariance 
methods  are  used.  (The  data  windowing  used  in  the  autocorrelation  method  spreads  some 
of  the  signal  power  into  the  noise  eigenvalues  making  it  difficult  to  separate  the  signal  and 
noise  subspaces  for  small  data  vectors).  However,  by  performing  the  eigen  decomposition 
on  the  estimated  autocorrelation  matrix  and  then  using  its  signal  subspace  equivalent  (i.e. 
the  smallest  N-M+l  eigenvalues  are  set  to  zero)  improved  signal  bearing  estimates  result. 

Since  the  deletion  of  the  noise  subspace  results  in  a  reduced  rank  autocorrelation 
matrix,  there  are  an  infinite  number  of  solutions  for  the  tap  weight  coefficient  vector.  For 
the  same  reasons  discussed  previously  in  section  3.0,  the  solution  vector  is  chosen  as  the 
one  with  the  minimum  norm.  The  improvement  in  the  DF  estimates  -  as  manifested  by 
the  suppression  in  spurious  peaks  -  using  these  techniques  and  shown  in  Figure  10  is 
obvious  when  compared  to  Figure  8. 
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FIGURE  11:  Bearing  estimation  accuracy  versus  model  order  for  the  case  of 

equi-amplitude  signals  (SNR  =  20  dB)  with  bearings  of:  (a)  40  and  100 
degrees,  and  (b)  40,  70,  100, 140,  and  160  degrees. 
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Figure  11  shows  the  effect  of  changing  model  order  on  the  accuracy  of  the  spectral 
estimator  for  the  identical  data  set  as  was  usea  in  Figures  3  and  9.  In  this  case  the  effects 
of  noise  have  been  significantly  reduced  which  has  the  effect  of  relaxing  the  restriction  on 
the  choice  of  the  optimum  model  order.  Kumaresan  and  Tufts  have  experimentally 
determined  that  a  model  order  of  p=(3/4)iVis  a  reasonable  choice  [8]  which  is  supported  by 
the  results  shown  in  Figure  11.  However,  in  using  the  modified  covariance  method,  a  choice 
of  p=(2/Z)N  might  be  more  suitable  since  this  results  in  the  maximum  rank  of  the 
autocorrelation  matrix,  and  correspondingly,  the  maximum  number  of  signals  that  may  be 
estimated,  without  any  significant  loss  in  resolution. 
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FIGURE  12:  Effect  of  the  choice  of  signal  number,  to,  on  the  estimated  DF  spectrum.  AR 
filter  order  p=10  (see  Figure  10(j)). 
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Up  to  this  point  it  has  been  assumed  that  the  number  of  signals  was  known  which 
would  allow  one  to  make  the  division  between  the  noise  and  signal  subspaces  in  equation 
5.2.  In  practice  this  is  an  unrealistic  assumption  and  some  method  is  required  to  determine 
this  value  since  spurious  estimates  can  result  when  the  incorrect  value  is  chosen  as  shown 
in  Figure  12.  One  method  that  has  been  proposed  is  to  examine  the  eigenvalues  to  try  and 
determine  the  division  between  the  larger  signal  subspace  eigenvalues  and  smaller  noise 
subspace  eigenvalues  (see  Figure  13).  This  is  analogous  to  examining  the  forward  prediction 
error  variance  to  determine  model  order  as  discussed  in  section  3.0,  and  the  problems  in 
doing  the  determination  this  way  are  identical. 
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FIGURE  13;  Eigenvalues  of  the  autocorrelation  matrix  used  in  Figure  12 


Using  a  similar  approach  as  discussed  in  section  3.0,  Kailath  and  Wax  have 
reformulated  the  AIC  and  MDL  algorithms  based  on  using  the  eigenvalues  of  the 
autocorrelation  matrix  [11].  The  results  are  given  by, 


and, 


AIC(m)  =  N  (p  -  m)  In 


+  m(2p-  m), 


(4.3) 


MDL(m)  =  N  (p  -  m)  In 


+  |m(2p-m)  ln(A) 


(4.4) 


The  optimum  signal  number  m  is  found  by  minimizing  the  above  functions.  The  values  A, 
are  the  eigenvalues  of  the  p*p  autocorrelation  matrix  which  is  formed  by  deleting  the  first 
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row  and  first  column  of  the  (p4T)*(p+l)  estimated  autocorrelation  matrix  (formed  using 
either  the  covariance  or  modified  covariance  method). 

Figure  14  shows  two  examples  of  how  both  the  AIC  and  MDL  algorithms  perform 
for  the  same  data  used  in  Figure  6  (section  2.0). 


Estimated  Signal  Number 


FIGURE  14:  Signal  number  determination  algorithms  for  a  (a)  two  signal  case,  and 

(bj  a  five  signal  case 


At  the  present  time,  little  is  available  in  the  open  literature  describing  the 
performance  of  these  algorithms.  To  provide  some  idea  of  their  performance  both  the  AIC 
and  MDL  algorithms  were  tested  against  sensor  data  simulated  for  a  16  element  DF  array 
over  100  trials  (the  number  of  signals  was  5).  The  method  of  generating  the  data  was 
identical  to  that  described  for  Figure  3(b)  and  Figure  7  in  section  2.0.  The  results  are 
shown  in  Figure  15.  In  this  particular  example  both  the  AIC  and  MDL  algorithms  were 
perfect. 


To  test  these  algorithms  against  a  more  difficult  signal  environment,  the 
simulations  were  re-run  with  signal  powers  randomly  chosen  to  be  between  10  and  40  dB 
greater  than  the  noise  power,  signal  phases  randomly  chosen  between  0  and  27t,  and  signal 
bearings  randomly  chosen  between  0  and  180  degrees.  Figure  16  (a)  shows  the  results. 
Under  these  conditons  the  performance  of  both  algorithms  was  severely  degraded  with  the 
correct  number  of  signals  (5)  being  determined  less  than  50%  of  the  time.  However,  in 
examining  the  cases  when  the  AIC  or  MDL  algorithm  failed  to  correctly  determine  the 
number  of  signals,  it  was  determined  that  in  most  of  these  cases  it  was  not  possible  to 
accurately  estimate  all  signal  bearings  from  the  DF  spectrum.  Consequently  the  trials  for 
Figure  16(a)  were  rerun,  except  trials  where  any  of  the  five  estimated  bearings  deviated 
from  the  true  bearing  by  more  than  5  degrees  were  discarded.  Figure  16(b)  shows  the 
results. 
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FIGURE  15:  Performance  of  the  model  order  determination  algorithms  using  the 

sensor  data  from  Figure  3(b),  i.e.  100  trials  of  a  16  element  system  in  a 
five  signal  environment. 
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FIGURE  16  Performance  of  the  model  order  determination  algorithms  as  in  Figure  15 
except  (a)  signal  powers,  phases,  and  bearings  were  randomly  chosen,  (b) 
same  as  (a)  except  trials  where  estimated  bearings  were  >  5  degrees  were 
ignored. 
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5.0  SUMMARY  AND  CONCLUSIONS 


In  this  report,  methods  to  determine  the  optimum  model  order  of  adaptive  filter 
techniques  (in  particular  autoregressive  filters)  applied  to  the  problem  of  radio  direction 
finding  were  evaluated.  In  filter  terms  model  order  is  the  number  of  filter  stages.  In 
direction  finding  (DF)  terms,  the  model  order  determines  the  antenna  subarray  size  for 
spatial  smoothing.  In  the  case  of  direction  finding  techniques  using  the  eigen  decomposition 
method,  the  related  problem  of  determination  of  the  number  of  signals  was  also  considered. 

Spatial  smoothing  is  a  technique  which  is  used  to  improve  direction  finding 
accuracy  in  the  presence  of  multipath  by  using  antenna  subarrays  (smaller  than  the  full 
array)  to  decorrelate  the  various  multipath  signal  components.  For  the  noncoherent  signal 
case  (uncorrelated  signals  transmitted  by  separate  sources),  this  technique  would  generally 
not  be  required  since  the  signals  can  be  decorrelated  in  time,  and  consequently  the  full 
array  is  used  (which  corresponds  to  a  model  order  of  N- 1  where  N  is  the  number  of 
antennas  in  the  DF  array). 

In  a  multipath  environment,  the  optimum  model  order  is  generally  not  known 
before  hand,  so  a  number  of  model  order  determination  methods  have  been  proposed  to 
tackle  this  problem.  In  this  report,  these  were  classified  according  to  the  autocorrelation 
matrix  estimator  used.  Three  classes  of  estimators  were  considered.  They  were:  (a)  the 
autocorrelation  method,  (b)  the  covariance  and  modified  covariance  method,  and  (c)  the 
eigen  decomposition  method. 

The  autocorrelation  method  is  not  appropriate  for  accurate  direction  finding 
purposes  since  adaptive  filter  DF  techniques  based  on  this  method  have  very  poor 
resolution.  However,  it  was  discussed  since  the  first  model  order  determination  algorithms 
were  proposed  for  this  method,  and  many  of  the  problems  encountered  are  common  to  all 
classes.  Four  algorithms  were  considered,  namely,  Forward  Prediction  Error  (FPE), 
Criterion  Autoregressive  Transfer  function  (CAT),  Akaike  Information  Criterion  (AIC), 
and  Minimum  Description  Length  (MDL).  The  performance  of  these  algorithms  is  quite 
poor,  typically  picking  model  orders  that  are  too  low. 

The  covariance  and  modified  covariance  methods  are  much  superior  for  direction 
finding  purposes  than  the  autocorrelation  method.  It  has  been  determined  both  analytically 
and  experimentally  that  a  model  order  of  between  N/ 3  and  N/2  is  suitable.  At  higher 
model  orders  spurious  peaks  begin  to  appear  in  the  DF  spectrum.  This  leads  to  a  severe 
degradation  in  the  bearing  accuracy  since  these  spurious  peaks  are  often  indistinguishable 
from  signal  peaks,  and  consequently  erroneous  bearing  estimates  are  made  using  these 
peaks. 


The  eigen  decomposition  method  enhances  the  signal  to  noise  ratio  of  the 
autocorrelation  matrix  estimated  using  either  the  covariance  or  modified  covariance 
method.  This  enhancement  suppresses  spurious  bearing  estimates  so  that  higher  model 
orders  can  be  used.  Values  of  between  2N/3  and  3N/4  have  been  proposed. 

One  difficulty  with  using  the  eigen  decomposition  method  is  that  the  number  of 
signals  is  required.  In  practice,  this  value  is  generally  unknown.  This  problem  is  closely 
related  to  the  model  order  determination  problem,  and  as  a  consequence,  both  the  AIC  and 
MDL  model  order  determination  algorithms  have  been  reformulated  for  this  problem.  In 
relative  terms,  the  performance  of  these  algorithms  is  much  superior  than  the 
corresponding  algorithms  tested  for  the  autocorrelation  method.  From  the  results  of 
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simulations  performed  for  this  report,  it  is  expected  that  the  performance  of  these 
algorithms  will  degrade  under  realistic  conditions  (i.e.  closely  spaced  signals  in  bearing,  low 
signal  to  noise  ratios,  etc.),  but  that  this  will  occur  at  the  same  time  as  the  accuracy  of 
eigen  based  DF  algorithms  also  degrade. 

In  general,  signal  determination  methods  used  in  conjunction  with  DF  methods, 
where  both  are  based  on  using  the  eigen  decomposition  method,  appear  to  offer  the  best 
performance  when  the  accuracy  of  the  DF  bearing  estimates  is  considered.  Model  order 
determination  methods  based  on  the  autocorrelation  method  were  found  to  perform  poorly 
and  are  considered  inappropriate  for  DF  purposes. 
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